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Properties of two algorithms for iterative reconstruction of SPECT images, LS-MLEM and LS-OSEM,are 
studied and compared with the ML-EM algorithm in this paper. By using projection data of heavy-noise, their 
effectiveness in improving SPECT image quality is evaluated. A phantom with hot and cold lesion is used in the 
investigation. The reconstructed images using LS-MLEM or LS-OSEM show that there is not a rapid increase in 
image noise,and the “best” estimate is assuming that the reconstructed images satisfy the statistical model. The 
major advantage of using LS-MLEM or LS-OSEM algorithm in SPECT imaging is in their ability to accurately 
control for heavy-noise. And LS-OSEM algorithm obviously improves the convergence rate. 
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I. INTRODUCTION 


In single photon emission computed tomography (SPECT) 
applications, the projection data of radiotracer activity are 
heavily contaminated by statistical noise. The noise-control 
methods of iterative image reconstruction and incorporat- 
ing statistical models can provide better noise properties 
than the conventional filtered back-projection method. The 
maximum-likelihood expectation-maximization(ML-EM) al- 
gorithm and ordered-subset variant (OSEM) algorithm [1, 2] 
are widely applied. However, the impact of statistical noise is 
a key problem in the image reconstructions, and it is crucial 
to estimate uncertainties of the reconstructed images [3]. 

In practical applications, firstly, the iterative algorithm is 
terminated before convergence to control the noise in re- 
constructed images, for higher-noise (low-count) specially. 
Only a few iterative algorithms with explicit multiplicative 
update equations can be applied. So many iterative al- 
gorithms, such as those based on principles of maximum- 
likelihood, penalized-likelihood and Bayesian maximum a 
posteriori(MAP) estimation, are applied to image reconstruc- 
tion. However, these iterative algorithms are nonlinear, and 
a theoretical formulation for noise propagation in image re- 
construction has been intractable. Whichever noise-control 
algorithm is used, at least one arbitrary parameter (often sev- 
eral) shall be chosen, but how to choose it (them) is a key 
issue. 

Secondly, those iterative algorithms are derived by assum- 
ing a low-noise (high-count) approximation, which means 
that the noise in the reconstructed images is much smaller 
than the mean images. Therefore, they become less accu- 
rate at higher-noise. As iteration number increases, noise 
variance in a reconstructed image tends to increase. So, a 
noise-controlling iterative algorithm in ML-EM reconstruc- 
tion stops before the ML point is reached [4, 5], and the final 
image obtained depends on the stopping point. 
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For statistical estimation problems, the best estimate can 
also be determined using the least squares (LS) and weighted 
least squares (WLS) criterion. The iterative algorithms pro- 
posed for LS or WLS image reconstruction include the op- 
timal step size (OSS) algorithms, the gradient descent (GD) 
algorithms, the conjugate gradient (CG) algorithms, etc. The 
primary difference among them is the way to determine the 
step direction. CG algorithms are used the most widely in 
SPECT image reconstruction [6-8]. The step direction de- 
pends on the previous step direction and the gradient direc- 
tions. The advantage of CG algorithms is that each new 
step dose not spoils the work of previous steps. However, 
this is not the case for LS/WLS-CG algorithms with a pre- 
conditioner, i.e. a transformation matrix. Performance of 
the LS/WLS-CG algorithms depends heavily on the pre- 
conditioner chosen. In general, the convergence rate of the 
WLS-CG algorithm is about ten times that of the ML-EM 
algorithm. Also, at large iteration numbers, the WLS-CG ex- 
hibits a faster increase in image noise than the ML-EM al- 
gorithm. The WLS-CG algorithm uses a Gaussian statistical 
model instead of the more correct Poisson model used in ML- 
EM algorithms. So, the WLS-CG algorithms are less accurate 
at higher-noise (low-count), and they are not easy to incorpo- 
rate the non-negativity constraint in the reconstruction, unlike 
ML-EM algorithms. According to standard least squares, the 
mathematical form of the LS criterion is an equation. Because 
of the equation elements contaminated by statistical noise and 
the large dimension of the equation, the close-form solutions 
of the equation are not often used in image reconstruction 
directly. The statistical noise is a Poisson model at higher- 
noise. Therefore, it is proposed that the equation based on 
standard least squares should be estimated using the ML-EM 
algorithms. Termed as LS-MLEM algorithm in this paper, it 
is advantageous in avoiding negativity and a faster increase 
in image noise at large iteration numbers. Mainly, the LS- 
MLEM based on least squares and ML-EM satisfies the Pois- 
son model and a low-noise approximation. 
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Il. LS-MLEM ALGORITHM 


In SPECT, the process of obtaining projection data is as 
follows: 


P=C), (1) 


where P is the vector with the elements representing photon 
counts measured in a given projection bin, Ais the vector with 
the elements representing photon counts emitted from a given 
voxel, and C is a system matrix with the elements c; repre- 
senting the relative contribution of J voxel to projection bin 
pi. The matrix C may model the effects of attenuation, de- 
tector response, and/or scatter, but these effects are ignored 
in this paper, so as to evaluate the result of reconstructed im- 
age based on LS-MLEM algorithm. Since the measurement 
of P is contaminated by noise, P is frequently modeled as 
a random process, and the elements are independent Poisson 
random variables. 
The mathematical form of the LS criterion is given by 


E(A) = ||P — CA|’. (2) 


The squared residual error F(A) shall be minimized. Accord- 
ing to standard least squares, we have 


C'7P = CC), (3) 


where C7 is transpose matrix of the system matrix C. 
Defining X = C'P, its elements x; are determined by 


I 
=> Gp FH 1,2) 4, (4) 
i=l 


where, J is the total number of given projection bin, and J is 
the total number of given voxel in a slice. P includes a Pois- 
son noise. Therefore x; in Eq. (4) is a linear superposition 
of the projection data with Poisson noise. Therefore, the vec- 
tor X, with its elements xj, also include Poisson noise. The 
close-form solution of Eq. (3) does not satisfy non-negativity 
constraint in the reconstruction. Signal-to-noise ratio SN Rx 
of the vector X is given by SNR, x I'/?SNRp, SNR, is 
signal-to-noise ratio of the vector P. Since the total number 
of projection bin is several thousands, the SNR, being ten 
times more than the SN R, can be approached. 

In Eq. (3), let H = C'C be a constant matrix based on 
system matrix, its elements hj; is determined by 


I 
hy = X cijca jl = 1a iJ. (5) 
t=] 


Because Eq. (3) is of lower-noise, ML-EM iterative algo- 
rithm based on statistical model is applied to estimate Eq. (3). 
Taking into account the Poisson statistics in the vector X, the 
LS-MLEM algorithm can be written as 


dj hijti 
Lj hig SD hua 


(6) 


k+l _ 
Aj = 
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Fig. 1. The true phantom. 
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Fig. 2. SN Rx versus SN Rp. 
As compared with Eq. (6), the conventional ML-EM algo- 
rithm of Eq. (1) is written as 


AK Cijpi 
etl = J why : (7) 
: Di Gij 3 D caf 


where k is the k® iteration. When the projection data includes 
higher-noise, the iteration for Eq. (6) is an iterating process 
with lower-noise, unlike Eq. (7) with higher-noise. It is well 
known that the ML-EM algorithm can get accurate image at 
lower-noise. So the noise propagated in reconstructed image, 
obtained using Eq. (6), can be controlled under lower-noise 
level. 


HI. COMPUTER SIMULATION 


To evaluate the properties of the LS-MLEM algorithms, 
we use a simulation phantom with hot and cold lesion. As 
shown in Fig. 1, the voxel is 3.00 mm and true image is of 
121 x 121 pixels, obtained using calculation of 2420 x 2420 
small pixels. The hot-to-background ratio of 5:lis assumed, 
and activity of the cold is zero. The radius of background is 
18.15cm, and the radius of hot and cold is 4.0cm. The dis- 
tance between hot and cold centers is 18.0cm. A low-energy 
(140 keV) general-purpose parallel-beam collimator (square 
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Fig. 3. Reconstructed images of the phantom in Fig. 1 (The LS-OSEM results will be discussed later). 
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Fig. 4. SNR versus iteration number (—, ML-EM; ---, LS-MLEM; ---, LS-OSEM). 


holes,of 2.24mm size, 4.90cm length, and 0.76 mm septa) 
is used. The SPECT head rotates around the phantom, with 
radius of rotation of 18.15cm. The effects of scatter, pene- 
tration and detector response function are ignored. Using the 
phantom, the simulated projection data are binned into 121- 
element arrays for each of the 120 views over 360°. Pois- 
son noise fluctuations are then added to the simulated projec- 
tion data. The projection data with SNR, = 5, 10, 20, 30, 
40, 50, 60, 70, 80, 90 and 100, i.e. the count/bin ratios of 


25, 100, 400, 900, 1600, 2500, 3600, 4900, 6400, 8100 and 
10 000 respectively, are used in evaluating the improved noise 
SN Rx. In Fig. 2, the SNR, data are compared with SN Rp 
data. The same projection data with noise are used to calcu- 
late the SN Ry and SN Rp. Fig. 2 indicates that the SN Rx is 
proportional to the SN Rp, but the SN Rx is almost 25 times 
larger than the SN Rp. The noise of vector X is remarkably 
improved, and the LS-MLEM is an iteration algorithm with 
lower-noise. 
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To emphasize LS-MLEM applied in projection data with 
higher-noise,projection data sets of SN Rp = 5, 10, 20 and 
30, are used in the reconstructions. The system matrix is ob- 
tained by calculating 2420 x 2420 small pixels. Both ML-EM 
and LS-MLEM iterative algorithms are used in reconstruct- 
ing the simulated projection data. Images of 121 x 121 pixels 
are reconstructed from 120 projections sampled over 360°. 
Results from the higher-noise data allow us to study proper- 
ties of the iterative reconstruction methods with LS-MLEM. 
An important property is its convergence rate. Also, the four 
higher-noise data sets allow comparison of LS-MLEM with 
ML-EM in noise improvement. And we have a better assess- 
ment of the improved reconstructed image quality. 


IV. RESULT AND DISCUSSION 


Figure 3 shows the reconstructed images, obtained using 
the ML-EM and LS-MLEM with Poisson noise added to 
the projection data, after the 20", 40", 90", 140", 500" and 
1000" iteration. 

Figure 4 shows the signal-to-noise ratio as function of iter- 
ation number for the ML-EM and the LS-MLEM algorithms, 
plotted by the solid line and dash dot line, respectively. The 
signal-to-noise ratio for the full, background and hot recon- 
struction field are respectively plotted from first to third col- 
umn. As Fig. 4 shows, for ML-EM, the maximum signal- 
to-noise ratio in full field is at 20, 40", 90" and 140" iter- 
ation at SNR, = 5, 10, 20 and 30 respectively; while for 
LS-MLEM, the signal-to-noise increases all the way to the 
end. 

Figure 3 shows that the noise artifacts in the reconstructed 
image, obtained using ML-EM, tend to increase with iteration 
number. This is more obvious as the SN Rp decrease which is 
consistent with Fig. 4. However, Figs. 3 and 4 show that noise 
in the LS-MLEM-reconstructed image tends to decrease with 
increasing iteration number where the artifacts do not change 
obviously. 

However, it is obvious that the convergence rate of the LS- 
MLEM is slower than the ML-EM. To increase the conver- 
gence rate, ordered-subset is applied to LS-MLEM algorithm 
(called LS-OSEM in this paper). The results, obtained us- 
ing LS-OSEM with 10 subsets, are plotted in Figs. 3 and 4. 
The convergence rate using LS-OSEM is improved obviously. 
Compared with ML-EM, the noise and artifacts from LS- 
OSEM do not change by much, though their increase can 
be seen in the reconstructed images with SNR, = 5 and 
10, at higher iteration number. A maximum signal-to-noise 
ratio is at 170" iteration for SNR, = 5. The result of 
SNR, = 20 and SNR, = 30 consist with LS-MLEM. As 
shown in Figs. 3 and 4, the SNR and quality of image, ob- 
tained using LS-OSEM, are better than the ML-EM, for hot 
field especially. 

To further evaluate the iterative algorithms in terms of im- 
age quality, Fig. 5 shows the line profile along axis of the cen- 
tral transverse cross section. The dotted line is of the phan- 
tom, while the other lines are defined the same as in Fig. 4. 
The line profile of ML-EM is at maximum signal-to-noise ra- 
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Fig. 5. Profile (---, true phantom, —, ML-EM; ---, LS-MLEM; ---, 
LS-OSEM). 
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tio, while the line profile of LS-MLEM and LS-OSEM are 
at the 1000" and 170" iterations, respectively. The horizon- 
tal axis is pixelnumber and vertical axis is activity. Note that 
Fig. 5 shows the image reconstructed by ML-EM algorithm 
higher noise than those by the LS-MLEM and LS-OSEM 
algorithm, and that the line profile of LS-MLEM and LS- 
OSEM algorithm is the “best” estimate of assuming that the 
line profile satisfy the statistical model. Figs. 3, 4 and 5 also 
indicate that the difference between ML-EM and LS-MLEM 
or LS-OSEM decreases with increasing SN Rp. 


V. CONCLUSION 


Using standard least squares, Eq. (1), of higher-noise, can 
be converted to Eq. (3), of lower-noise. The SNR, is ten 
times more than the SN Ry. This can provide an algorithm of 
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